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Abstract. We discuss the problem of heat conduction in quantum spin chain mod- 
els. To investigate this problem it is necessary to consider the finite open system 
connected to heat baths. We describe two different procedures to couple the system 
with the reservoirs: a model of stochastic heat baths and the quantum trajecto- 
ries solution of the quantum master equation. The stochastic heat bath procedure 
operates on the pure wave function of the isolated system, so that it is locally and 
periodically collapsed to a quantum state consistent with a boundary nonequilib- 
rium state. In contrast, the quantum trajectories procedure evaluates ensemble 
averages in terms of the reduced density matrix operator of the system. We apply 
these procedures to different models of quantum spin chains and numerically show 
their applicability to study the heat flow. 

1 Introduction 

The derivation of Fourier's law of heat conduction from the microscopic dynamics, without 
any ad hoc statistical assumption, is one of the great challenges of nonequilibrium statistical 
mechanics [1]. This problem, in spite of having a long history, is not completely settled: Given a 
particular classical, many-body Hamiltonian system, neither phenomenological nor fundamental 
transport theory can predict whether or not this specific Hamiltonian system yields an energy 
transport governed by the Fourier law J = -kVT, relating the macroscopic heat flux to the 
temperature gradient VT [2]. 

At the classical level a great amount of work has focused on the relation between the chaotic 
character of the microscopic dynamics and the normal macroscopic transport [3-18], (see also 
Ref. [19] for a recent review ). The general picture that emerges from these investigations is 
that deterministic chaos appears to be an essential ingredient required by transport theory. 
However, strictly speaking, the exponential instability that characterizes the chaotic dynamics 
is neither sufficient [7] nor necessary [10, 16] for the validity of Fourier law. What has been 
shown in Ref. [16] is that the diffusive behavior, which is at the root of normal heat transport, 
can be obtained even if the Lyapunov exponents are zero. This constitutes a strong suggestion 
that normal heat conduction can take place even without the strong requirement of exponential 
instability. 

At the quantum level, the question whether normal transport may arise from the underlying 
quantum dynamics is a controversial issue [20-26]. This is mostly because it is not clear how 
to describe the transport of energy or heat from a microscopic point of view. In analogy to 
classical systems, a quantum derivation of the Fourier's law calls directly in question the issue 
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Fig. 1. Schematic representation of a finite one-dimensional quantum spin chain, coupled to external 
heat reservoirs at different temperatures. 



of quantum chaos [37]. However, a main feature of quantum motion is the lack of exponential 
dynamical instability [27]. This fact may render very questionable the possibility to derive the 
Fourier law of heat conduction in quantum mechanics. Thus it is interesting to inquire if, and 
under what conditions, Fourier law emerges from the laws of quantum mechanics (for a recent 
review of the microscopic foundations of the quantum Fourier's law see [28]). 

To investigate the problem of quantum heat transport one has to deal with a finite open 
system connected to heat baths. This fact renders the problem to obtain a derivation of the 
quantum Fourier's law extraordinary difficult. In this paper we present two complementary 
methods to study the problem of heat conduction in quantum spin chain models. We consider 
one-dimensional quantum spin chain models like the one schematically shown in Fig. [TJ Each 
spin in the chain interacts with its neighbours and possibly with an external field. The range 
and type of interaction does not limit the applicability of our methods. Furthermore, the spin 
chain is in contact with two heat reservoirs. The Hamiltonian of the system can be written as 

H system ~t~ HreservoirsiT^L^^R^j ~\~ Hint j (1) 

where Hint represents the coupling between the system and the heat reservoirs. 

In recent years the transport of heat in quantum spin chains has been the subject of intense 
experimental [31] and theoretical investigations [20-26,32-36]. In analogy to classical systems 
it has been found that non integrable quantum spin chains show normal heat transport, while 
integrable chains lead to ballistic transport and thus, to divergent heat conductivity. In par- 
ticular, in Ref. [25] the validity of the quantum Fourier's law has been linked to the onset of 
quantum chaos, in the statistical sense of Random Matrix Theory. Moreover, in Ref. [25] an 
intermediate case for which the spin chain is neither integrable nor chaotic was shown to have 
divergent heat conductivity. Up to our knowledge, this is the only example for which a quantum 
spin chain with intermediate statistics leads to an abnormal heat transport (see also [26] where 
almost integrable models have been studied). In spite of the general evidence, in Ref. [22] the 
Fourier's law was numerically confirmed for a class of integrable, albeit small, quantum spin 
chains. 

In these investigations, particularly those that are numerical, the thermodynamical limit at 
which the heat conductivity is formally defined is difficult, if not impossible to address. This 
limitation is emphasized when one deals with non integrable systems for which a theoretical 
perturbative analysis is not possible. So far the most popular theoretical framework to study 
heat transport is the Green-Kubo formula [29] , derived on the basis of linear response theory. 
However, the use of GK formula for heat transport requires ad-hoc statistical assumptions that 
lack of a microscopic justification of its applicability (see e.g., [28] and also [30] for a plausible 
derivation of a Green-Kubo formula for heat transport). A second standard treatment is based 
on the Quantum Master Equation (QME). However, this method involves the calculation of 
the reduced density matrix of the system, thus limiting numerical investigations to relatively 
small system sizes. 

Given this state of affairs, the development of numerical methods that can deal with both 
integrable and non integrable spin models and that are amenable to study large system sizes is 
very desirable. In the present paper we present two methods recently introduced to study the 
transport of heat in quantum systems. In section [21 we discuss a model of stochastic quantum 
heat baths where Hi nt consists of a stochastic time-periodic local perturbation [25]. This per- 
turbative term acts on the, otherwise, pure state of the system. We show that this procedure 
leads to a well defined local temperature and discuss the validity of Fourier's law in a Quantum 
spin-1/2 Ising chain in a tilted external magnetic field. This method allows to study finite but 
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large spin chains as its implementation only requires the knowledge of the pure vector state of 
the system. In Section [3] we present a QME in Lindblad form that can appropriately describe 
nonequilibrium steady states [36]. Furthermore, we use the Monte Carlo wave function formal- 
ism to study the heat transport in the spin-1/2 Heisenberg chain. In Section |4] a discussion of 
the results and applicability of our integration methods is presented, followed in Section O by 
our final remarks. 



2 Stochastic Quantum Heat Baths 

In this Section we describe the implementation of a stochastic model of quantum heat baths. 

We consider a one-dimensional quantum spin-1/2 chain of length N. The Hamiltonian of 
the system can be written in general as in Eq. [1] Each spin in the chain is coupled to its 
neighbours and possibly to an external field. Furthermore, let us consider that the leftmost 
(so) and rightmost (sat_i) spins are coupled via the Hi nt term, with external ideal heat baths 
at temperatures Tl and Tr respectively. A schematic representation of this general model is 
shown in Fig. [TJ 

The aim of this procedure is to focus on the evolution of the state of the spin chain, avoiding 
to introduce any particular model for the heat reservoirs. Therefore, we only assume that the 
heat baths act locally on the state of the respective spin so that the state of so and of sn-i are 
thermal states at the respective temperatures. This model for the reservoirs is analogous to the 
stochastic thermal reservoirs used in classical simulations [15] and we thus call it a stochastic 
quantum heat baths. 

To be precise, we work in the representation basis of o~ z . Furthermore, we use units in which 
Planck and Boltzmann constants are set to unity h = fee = 1 In this the wave function at time 
t can be written as 

W{t))= C SO,Sl,-,SN-l( t )\SO,Sl...S N - 1 ) , (2) 

So,5l,...,SJV_l 

where s n — 0, 1 represents the up, down state of the n-th spin, respectively. The wave function 
at time t is obtained from the unitary evolution operator U(f) = exp(— iTLt). The interaction 
with the reservoir is not included in the unitary evolution. Instead, we assume that the spin 
chain and the reservoir interact only at discrete times with period r at which the states of the 
spins so an d sn-i are stochastically reset. The evolution of the wave function from time t to 
time t + t can be represented as 

h/>(i + r))= S(T L ,T R )V( r )m)) , (3) 

where S(Tl,Tr) represents the unitary stochastic action of the interaction with the left and 
right reservoirs at temperatures Tl and T# respectively. 
The action of S(Tl, Tr) takes place in two steps: 

(i) A local measurement of the state of the spins coupled to the heat baths is performed. Then 
their state collapses to a state (sq 1 s* l _ 1 ) with probability 

P(4' S *N-l)= \ C s* ,°U..,SN-2,s' N J 2 ■ ( 4 ) 

Si,...,SiV-2 

Numerically this means that we put all coefficients C SOjS1i ... iSjv _ 1 with (so,sat_i) =/= { s qi s *n-i) 
to zero. Afterwards, the wave function is renormalized. 

(ii) The new state of the edge spins is stochastically chosen: sq and sn^i are set to down, (up) 
state with probability /i,(l — fi). The probability fj.([3) depends on the canonical temperature 
of each of the thermal reservoirs: 
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where (3 = 1/T is the inverse temperature. This simulates the thermal interaction with the 
reservoirs. 

This interaction thus (periodically) resets the value of the local energy of the spins in contact 
with the reservoirs. This information is then transmitted to the rest of the system during its 
dynamical evolution and relaxation towards equilibrium. Therefore, the value of t controls the 
strength of the coupling to the bath. We have found that, in our units, r = 1 provides an 
optimal choice. On one hand, large values of r correspond to a "microcanonical" situation: if r 
is much larger than the relaxation time of the system then the spin chain behave as an isolated 
system, leading to a equidistributed mean energy. On the other hand, a strong coupling with the 
reservoirs, i.e., t <C 1, the state of the system freezes due to the quasi-continuous monitoring, 
a situation similar to the Zeno effect. We have checked that intermediate values of r lead to 
qualitatively similar results. 

Also, note that our method does not correspond to a stochastic unraveling of a QME. 
Here the evolution of the system between two consecutive interactions with the baths is purely 
unitary and not dissipative as for the QME in Section \3\ 

The use of stochastic quantum heat baths has the advantage that it only requires the 
calculation of the vector state of the system (J3)). This allows to computed time averages for 
spin chains longer than the sizes that can be studied with other methods. 

Finally, we remark that the stochastic quantum heat baths method does not depend on the 
range and type of the interaction. 

2.1 Fourier's law in the Ising chain in a tilted magnetic field 

In this Section we discuss the heat transport in a Ising chain of N spins 1/2 with coupling 
constant Q subject to a uniform magnetic field h = (h x ,0,h z ), with open boundaries. The 
Hamiltonian reads 

L-2 , 

H= £fln + -(<7 1 +<7 r ) • (6) 

n=0 

where H n are local energy density operators 

H n = -Qa*a* +1 + — ■ (a n + cr n+1 ) , (7) 

and o\ = h • <To/h, a T = h • a^-i/h are the spin operators along the direction of the magnetic 
field of so and sn-i respectively. The operators a n = (a^,a^,a^) are the Pauli matrices for 
the n-th spin, n = 0, 1, ... TV — 1. The direction of the magnetic field affects the qualitative 
behavior of the system: it is integrable for h z = and not integrable otherwise. When h z is 
of the same order of h x quantum chaos sets in [25] . Using the stochastic quantum heat bath 
formalism, in Ref. [25] it was shown that the validity of the Fourier's law is related to the 
onset of quantum chaos. In what follows, we consider a chaotic chain (with parameters Q = 2, 
h x — 3.375, h y — and h z = 2) and discuss the establishment of local thermal equilibrium in 
terms of the density matrix operator. Furthermore, we also recall some of the results presented 
in Ref. [25] concerning the validity of the Fourier's law. 

In order to apply the stochastic formalism to this model one needs first to rotate the 
state of the edge spins to the direction of the external field h, so that the local measure- 
ment described above (z) is meaningful. This is done by rotating the wave function by the angle 
a = ta,TT 1 {h x /h z ) to the eigenbasis of the components o\ and <r r , that is \tp) — > e~ la ^ cr ° +cri '- 1 " 2 \tp) 
After the stochastic reset of the edge spins the wave function is rotated back to the u z basis, 

To integrate the unitary evolution U(t)\t/j) of the system we have used an accurate high 
order split-step factorization of the unitary evolution operator described in Ref. [38] . We then 
consider different random quantum trajectories of the randomly chosen initial wave function 
\ip(0)) of the system. The state is then evolved for some relaxation time r ro i after which it is 
assumed to fluctuate around a unique steady state. Measurements are then performed as time 
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Fig. 2. Effective temperature at the bulk /3 _1 as a function of the nominal value of the temperature 
of the baths T. (3 was obtained from a best fit to a exponential of p n (E n ) in the central symmetry band 
for the chaotic chain of length N = 7. The dashed line corresponds to the identity. 



averages of the expectation value of suitable observables. We further average these quantities 
over the ensemble of "quantum trajectories" . 

In order to test the effectiveness of the coupling between the stochastic quantum baths and 
the system, we have computed the time averaged density matrix of the system 

p= ]im f |V(s)><^(s)lds , (8) 

where ip( s ) is the state of the system at time s. Quantum statistical mechanics postulates that 
if a system, described by a Hamiltonian H , is in a thermal equilibrium state at temperature 
then in the energy basis (H\<p n ) — E n \<p n )), the density matrix operator is 

e -/3B„ 

(<frn\p\<f>m) = — "g — $m,n , (9) 

where Z — ^2 n e~ f3En is the canonical partition function. 

We have performed equilibrium simulations, i.e., Tl = Tr = T and computed the time 
averaged density matrix. We have found that p is diagonal in within numerical accuracy. More- 
over, p n = (cf) n \p\(j> n ) is an exponential function of E n inside each symmetry band. This verifies 
Eq. [9] and thus, that the system reach a canonical equilibrium. Furthermore, from a best fit 
to a exponential of p n (E n ) one can extract a value for the effective temperature at the bulk 
of the system. In Fig. [2] we show the obtained effective bulk temperature /3 _1 as a function of 
the nominal value of the bath's temperature T. At low temperatures the effective temperature 
at the bulk saturates to a constant which, together with the energy profile E n , is in principle 
determined by the ground state. However, we have found that the saturation value does not 
correspond to neither the ground state of the two-body energy density operator H n nor to 
the ground state of the many-body operator Tt. At high temperatures (T > 5), the spin chain 
thcrmalizes to exactly the same temperature as that set by the stochastic heat baths. This 
strongly supports the effectiveness of our bath model. Moreover, this also suggests the validity 



6 



Will be inserted by the editor 



25 



20 



y 15 



10 



Opr 



'I I 1 



" 3 0.2 0.4 0.6 0.8 1- 

n/N 



+ + + ++ ++ + + 



J_ 

20 



J_ 
18 



J_ 
16 



J_ 
14 

N~ 



J_ 
12 



J_ 
10 



Fig. 3. Size dependence of the energy current for the chaotic chain with Tl = 5 and Tr = 50. The 
dashed line corresponds to a 1/N scaling. In the inset, energy profile computed from the time average 
of the expectation value of the energy density operator E n = (H n ). 



of the ergodic property, namely the averages over the canonical ensemble are equivalent to time 
averages. 

Out of equilibrium, we have verified that at high temperatures, the local temperature ob- 
tained from time average of the reduced density matrix of subsets of few spins centered around 
the n-th spin coincide with time averages of the expectation value of the energy density op- 
erator (H n ). However, we have observed that out of equilibrium the local temperature of the 
edge spins may not coincide with the nominal value of the temperature of the corresponding 
bath. This energy jumps are commonly observed in quantum and classical systems. They can 
be understood as the result of an thermal resistance of the particular contact model. 

Following Ref. [25], one can define a local energy current operator for the spin model of 
Eq. [H] as 

J n = h x Q «_! - ct*, l<n<iV-2, (10) 



that is consistent with the conservation of energy at the bulk. In Fig. [3] we show the heat 
conductivity k — J/VT as a function of the size N of the chaotic chain for sizes up to N = 20. 
The mean current J was calculated as an average of ( J n ) over time and over the N — 8 central 
spins. Three spins near each bath have been discarded in order to be in the bulk regime. For the 
particular choice of the energy density operator |(7j), its averaged expectation value is related to 
the local temperature as (H n ) oc — 1/T [25]. The temperature difference was thus obtained as 
AT = —1/(Hn-5) + 1/(1/3). For large iV the heat conductivity of the chaotic chain converges 
to a constant value, thus confirming the validity of the Fourier's law. In [25] it was also shown 
that for the integrable and an intermediate chains the Fourier's law does not hold as k diverges. 
The results of Ref. [25] supports the relation between a normal transport and the onset of 
quantum chaos. 
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3 Quantum Master Equation for Nonequilibrium States 

In this Section we present an alternative way of modeling a nonequilibrium transport scenario. 
We consider a spin-^ chain of length N with homogeneous and isotropic nearest-neighbor 
Heisenberg interaction V and a contribution H cyit due to the interaction with an external 
uniform magnetic field in z-direction. The Hamiltonian of the spin chain reads 

H = # cxt + V , (11) 

N-l N-l n 

H^=J2 H n =Y,J<> ( 12 ) 

N-2 N-2 
V = 2J V n ,n+1 = A 2J a n ' CT n+ l . (13) 

The energy contribution due to the external field of strength J? is assumed to be an approxi- 
mately conserved quantity. This is well justified if Q is large compared to the coupling strength 
A. By making use of a discretized version of the continuity equation, the following form of the 
energy current operator is found 

Jn,n+1 = i[V n ,n+l, H„] . (14) 

The central ideas of the following considerations stem from the theory of open quantum 
system, commonly used to study the properties of systems in thermal equilibrium, or to account 
for environmental effects in otherwise perfectly coherent quantum dynamics. Our aim is to 
describe the nonequilibrium steady state, that arises when the chain of spins is in contact with 
two thermal baths at different temperatures, leading to a steady energy current that will flow 
from the hotter towards the colder heat bath through the chain of spins. 

Our method is based on a Markovian quantum master equation for the reduced density 
operator ps of the spin chain, which is designed to model energy transport under a thermal 
gradient. In addition to the usual assumption of weak system-bath interactions, and to neglect 
all memory effects (Born/Markov approximation), we require the internal interactions between 
spins to be weak. By further assuming not too low temperatures of the baths the following 
QME is well justified (see [36], and references therein): 

±ps(t) = -i[H, Ps (t)}+V h (p s (t))+V K (p s (t)) (15) 

2 - 

^l(ps(*)) - £ (tl)m (F kP8 F} - -\F}F k , p s } + ) (16) 
k,i=i 

F x = a+, F 2 ee a - . (17) 

The dissipator of the right heat bath 2?r is defined correspondingly. The coefficient matrices 
7l/r read 



where /3l/r refers to the reciprocal temperature of left /right bath respectively, Ar controls 
the system-bath coupling strength and I(f2) denotes the spectral density of the bath, that we 
choose Ohmic. A remarkable property of Equation (fT5| is that it can be brought into Lindblad 
form [40, 41] 

-ps(t) = -ilH, Ps (t)}+J2a k (L k ps4 - i[4L fc ,ps]+) (19) 



by diagonalizing the coefficient matrices 7l/r • Therefore Equation (|15p can be treated with 
the Monte Carlo wave function technique that will be briefly described in the following Section. 
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3.1 Monte Carlo wave function method 

As the solution of QME's is rarely available analytically one is most likely confronted with 
the numerical time propagation of the reduced density operator. Standard methods for the 
numerical solution of linear differential equations like Runge Kutta solvers can readily be applied 
for small systems, but fail if the dimension of the Hilbert space d becomes large in view of 
computer memory limitations. 

Since the early 1990's the Monte Carlo Wave Function (MCWF) technique (also known as 
the quantum jump approach), has become increasingly popular. Introduced in the context of 
quantum optics [44, 45] the MCWF technique has also been used in quantum state diffusion 
unraveling of the QME [46] and in the study of the stability of quantum algorithms [47,48]. 
The basic idea is to depart from a statistical treatment by means of density operators and 
turn to a description in terms of stochastic wave functions. One might be tempted to think of 
single systems being continuously monitored, but this analogy is sometimes criticized. Never- 
theless, the MCWF technique has proved to be a powerful method for the numerical solution 
of Lindblad-type QME's. 

The Lindblad QME (fT9|) can equivalently be formulated in terms of a stochastic Schrodinger 
equation (SSE) 

■h 



|#) = (# efI + |) |</>) dt+Y, f ~|= - 1 j IV-) dn k (20) 



deterministic evolution 



jump-like evolution 



describing a piecewise deterministic process in Hilbert space [39] , a solution of which is is called 
a realization. The first term in Eq. (|20p corresponds to a deterministic time-evolution due to 
an effective non-Hermitian Hamiltonian given by 

HeS = -m-^J2 L t L >° ( 21 ) 
fe=i 

whereas the second term in Eq. (|2H)) refers to jump-like, stochastic evolution induced by the 
jump operators 

Jk = \fuk~L k . (22) 

The Poisson increments dn k € {0, 1} obey the following statistical properties 

E(dn fe ) = Pk dt , (23) 
dn k dni = 5 M dn k . (24) 

E(-) stands for the expectation value, whereas p k denotes a jump rate given by 

Pk = || Jk N>> f , (25) 

and is therefore time-dependent, p = J2 k Pk refers to the total jump rate. The SSE (|20p is an 
equivalent formulation of the Lindblad QME (fT9|) insofar as 



E(Mt)> = PS® , (26) 

given that E(|V>(to)) (^(to)|) = ps(to)- Thus the expectation value of an observable A at time t 
can be estimated through 

-. m 

(A) (t) = Tr{A Ps (t)} ~ - £ (Mt)\A\Mt)) (27) 



m 
fe=l 



in a finite ensemble of m realizations of (|20p to arbitrary precision. This is of huge practical 
importance, as one deals with wave functions with O(d) elements instead of density operators 
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having 0(d 2 ) elements, where d is the dimension of the Hilbert space under consideration. 
Furthermore, if one is interested in the stationary state, ensemble averages can be replaced 
by time averages and one single realization suffices to determine stationary expectation values 
[42,43] 

1 T 

(AY™=T±{A$"}~ ^-^J2(i,(t k )\Am k )), t k ^t + kAt. (28) 

The MCWF approximation to the exact solution depends on the number of quantum tra- 
jectories that are considered, as well as the integration time of each trajectory. We have found 
that one single trajectory, integrated during a sufficiently long time, results in a very good 
approximation. Nevertheless, following [43], we consider i = l...m different realizations, so 
that, an estimate of Tr{Ap| tat } can be obtained as the sample mean of the (j4)f at 's. Moreover, 
considering several quantum trajectories enables us to get a measure for the statistical error 
5(A) Btat /y/m, where S(A) stat is the standard deviation of the (A)f tat 's . The results presented 
below were obtained from m = 4 different realizations, integrated for long times of the order of 
fO 5 to 10 7 . 

We now describe the procedure to obtain a realization of the SSE l|20[) . 
Starting from an initial state IVK^o)) we proceed as follows: 

1. Draw a uniformly distributed random number r € [0, f]. 

2. Perform the deterministic time evolution | "4>(t)) — e^ - * '* 1 " \ip(to)) until t = tj, which is 
determined by ||V>(ij)ll 2 = r i f° r some r < 1. 

3. Normalize the wave function \ijj) — > | ijj) / 

4. Choose randomly a particular jump k with respective weight p k . 

5. Carry out the map \ip) — > J k / \fpk~- 

6. Set to — > t and return to step f . 

This procedure terminates when t = fg n , where ig n is the desired final time. 

The jumps occur at random instants of time tj, which are determined through HV'fe)!! 2 = r i 
by virtue of the second step of the simulation procedure. Assuming a uniform time discretization 
At in a simulation, jumps happen only at multiples of At, which causes an error of order O(At). 
Therefore At has to be chosen with care. For a more detailed account on the MCWF method, 
the reader is referred to [39] . 



3.2 Fourier's Law in the Heisenberg chain 

We consider the Heisenberg chain of spin-^ with Zeeman contribution, described by the Hamil- 
tonian (|11|> . as an example of an integrablc chain. As we have mentioned, integrability is com- 
monly associated with diverging transport coefficients and hence ballistic transport behavior. 
However, in Ref. [22], normal transport was observed for an integrable, albeit small, spin chain. 
In Fig. [4] we show the energy profile along the Heisenberg chain in the stationary state. A clean 
finite temperature gradient is observed. 

We explicitly compute transport coefficients here for Heisenberg chains of up to 12 spins, 
by making use of the MCWF technique, taking time averages over single trajectories. We focus 
on the scaling behavior of the heat conductivity n with the size of the chain N . We define k as 
the ratio of two measures, namely the stationary energy current within the chain, in terms of 
the current operator given in Eq. 114[) and the stationary difference in the local energy of the 
innermost pair of spins 



(Ji, 2 ) stat 



(29) 



AT^^iH^-Hv}^, v ee ^ - mod(A, 2) . (30) 

In the upper panel of Fig. [5] it is shown that the stationary current depends linearly on the 
reciprocal chain length N~ x . The same result has been observed earlier in Ref. [22] on the 
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Fig. 4. Energy profile in a Heisenberg chain of N = 10 spins- 1 . A finite gradient is associated with 
diffusive transport behaviour. The data points refer to time averages of single stochastic realizations 
up to a final time ig n = 10 5 in units of J? -1 . System parameters: (/3l = 0.41, /3r = 1.39, A = Ab = 
0.01, Q = 1) . 



basis of a similar QME, but for chain lengths of N < 6. Extrapolation of a linear best fit of 
the data points in the upper panel of Fig. [5] strongly suggests that, in the limit of infinite N 
the stationary energy current remains finite. In contrast to the diffusive behaviour observed for 
small systems, this is expected to wane at the thermodynamic limit N — ► oo . 

In a real material the heat conductivity n is a bulk property. When transport is normal, 
k converges towards a constant value with an increasing size of the system, once finite size 
effects are negligible. However, our results show no sign of a convergence for k in the range of 
our computational abilities. Linear extrapolation of (Ji t 2) stat and AT v _\ tV rather predicts a 
divergence of the so defined heat conductivity with N. This agrees with a ballistic transport in 
the integrable systems. 



4 Discussion 



In the previous sections we have described two formalisms to study the dynamics of quantum 
spin chain models in thermal nonequilibrium states. In this section we discuss their limitations 
and applicability . 

Both methods successfully generate a nonequilibrium state in the bulk of spin chain models. 
In the stochastic baths method the information of the model enters in the precise definition of 
the local measurement periodically performed on the spins in contact with the baths. For the 
QME, the particulars of the system are defined in the decay rates of Eq. fT8l 

However, the two methods have a very different character. The derivation presented in 
Section [3] obtains a QME of Lindblad type, appropriate to describe a quantum state for which 
the temperature field is not uniform. The fact that Eq.[l5]can be written in Lindblad form makes 
possible to compute averages of thermodynamical observables using the Monte Carlo wave 
function formalism, in which each quantum trajectory corresponds to a stochastic unraveling 
of the QME (15]). 

The model of Section [2] positively neglects any particular model for the heat baths and 
defines a procedure by which the system stochastically dissipates at its boundaries. As so, this 
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Fig. 5. Scaling behavior of heat conductivity. The different panels show the energy current (top), the 
temperature gradient (middle) and heat conductivity k (bottom) in the stationary state as a function 
of the reciprocal chain length iV -1 . The results where obtained by the MCWF method. In the upper 
and middle panels the dashed line corresponds to a linear best fit. In the middle panel only data 
points with N > 8 were considered. The error bars in k where obtained from error propagation in 
linear approximation. System parameters: (5l = 0.25, (3r — 0.5, A = Ab = 0.01, S2 = 1 . Simulation 
parameters: At = 1, t Rn = 10 7 (N > 9), t Rn = 10 8 (iV < 8) . 



formalism does not corresponds to a stochastic unraveling of a master equation for the density 
matrix operator. The evolution of the system between two consecutive interactions with the 
baths is purely unitary and not continuously dissipative as in the QME. The stochastic method 
is an approximation of the QME. 

The use of the MCWF technique to obtain averaged expectation values make possible to 
numerically study larger systems than with other methods. Without this, the solution of Eq.fTBl 
involves the integration in full Liouville space of dimension 2 2N , limiting numerical investiga- 
tions to small system sizes (typically of N < 6). The price to pay is that the Lindblad QME 
(fl~5|) . is valid only for weakly coupled spins. This limitation is particularly relevant for chaotic 
systems as typically, chaotic behaviour is exhibited above a critical interaction strength. Despite 
its limitations, Eq. [15] is, up to our knowledge, the only rigurous QME of Lindblad form that 
is adapted to study systems in a nonequilibrium state. 

The stochastic baths method does not present this limitation. The strength of the spin inter- 
action determines the internal relaxation time of the system. Since the frequency r _1 at which 
the system dissipates in the heat baths is a free parameter, one can always find an appropriate 
value of t for which, the stationary nonequilibrium state is established. Moreover, the stochas- 
tic baths can be generalized to consider more general situations, like e.g. the coupling with 
thermo-magnetic baths, for which thermodynamic cross effects can be studied. Nevertheless, 
the lack of a model for the baths limits a precise interpretation of the physical dissipation. 



5 Conclusions 

We have presented two complementary methods to study heat transport in quantum spin 
chains. The first is based on a stochastic procedure by which, the state of the subsystems that 
are coupled to ideal heat baths is consistent with a global nonequilibrium state. The second is 
based on a QME in Lindblad form. 
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The stochastic baths method does only require the integration of the pure wave function of 
the system. The Lindblad QME (I15|) can be integrated by means of Monte Carlo wave function 
techniques. Therefore, the use of any of these methods allows us to study longer spin chains 
than with other methods. This is particularly relevant to the study of nonequilibrium states 
as the quantities that determine the transport properties are formally defined in the limit of 
infinite volume. 

We have shown the application of these two methods to study the validity of the quantum 
Fourier's law in a chaotic and an integrable quantum spin chain. As generally observed, we have 
obtained that the Fourier's law holds for the chaotic chain, while for the integrable chain, the 
heat conductivity diverges with the size of the chain. 

It would be interesting to compare the nonequilibrium state generated by the two methods 
presented here. An investigation in this direction will appear elsewhere [49]. 

C.M.-M. acknowledge a Lagrange fellowship from the Institute for Scientific Interchange Foundation. 
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